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ABSTRACT 

The  solution  of  the  sparse  algebraic  system  for  a 
point  controlled  elliptic  distributed  parameter 
system  by  a multiple  shooting  and  sweep  algorithm 
Is  discussed.  The  multiple  shooting  and  sweep 
algorithm  enhances  the  convergence  rate  of  the 
ordinary  shooting  method  while  achieving  a sig- 
nificant reduction  In  the  dimension  of  the  linear 
equations  to  be  solved. 

The  optimization  algorithm  for  a special  problem  In- 
volving the  minimum  cost  selection  of  source  In- 
tensities with  the  state  satlsfing  a specified  con- 
straint set  Is  presented.  An  application  of  the 
techniques  Involving  the  analysis  and  management 
design  for  water  quality  control  In  Corpus  Christl 
Bay,  an  estuary  on  the  Texas  Gulf  Coast  In  the 
United  States,  Is  discussed.  The  resulting 
algorithm  is  more  than  twice  as  fast  as  a corre- 
sponding successive  overrelaxation  method. 

INTRODUCTION 

The  control  of  a class  of  distributed  parameter 
systems  which  are  described  by  elliptic  partial 
differential  equations  specified  over  a two  dimen- 
sional spatial  domain  and  satisfying  mixed  Dlrlchlet 
and  Neumann  type  boundary  conditions  (1)  Is  con- 
sidered. The  controls  for  the  system  are  defined  by 
parameters  whose  Influence  Is  exerted  polntwlse. 
Systems  of  this  type  are  encountered  for  example  In 
estuarine  water  quality  problems  (2)  where  the 
point  controls  correspond  to  pollution  discharges 
Into  estuaries. 

A method  of  solving  the  finite  difference  approxima- 
tion for  the  boundary  value  problem  baaed  on  a 
multiple  shooting  (3)  technique  is  presented.  It  Is 
equivalent  to  a Newton's  solution  of  the  system 
equations  generated  when  the  boundary  value 
problems  are  re-expressed  as  Initial  value  problems. 
The  method,  which  converges  for  the  elliptic  system, 
has  been  demonstrated  to  achieve  convergence  In 
two  example  problems  after  less  than  10  Iterations 
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and  Is  more  than  twice  as  fast  as  the  corresponding 
successive  over-relaxation  (SOR)  solution  methods 
(4,5).  The  order  of  linear  algebraic  equations  to 
be  solved  Is  the  maximum  of  N or  M where,  N and 
M are  the  number  of  grid  points  defined  In  the 
coordinate  directions  for  the  approximate  problem. 
This  compares  favorably  with  the  N X M order 
system  usually  required  In  the  other  methods. 

The  determination  of  optimal  cost  removal  policies 
for  discharges  by  nonlinear  programming  Is  discussed 
when  the  spatial  concentration  defined  by  the  system 
equations  are  required  to  satisfy  a set  of  state 
variable  Inequality  constraints.  The  sensitivity  for 
the  system  state  to  the  parametric  controls  Is  shown 
to  satisfy  a monotonlclty  condition.  This  mono- 
tonlclty  may  be  combined  with  convexity  of  the  cost 
functions  to  reduce  the  order  of  the  resulting  non- 
linear programming  problem  by  solving  a set  of 
problems  defined  for  relaxed  consistent  state  variable 
constraints  that  converge  to  the  desired  constraints. 

The  solution  technique  Is  applied  to  a study  of  the 
spatial  distribution  of  phosphorus  in  Corpus  Christl 
Bay,  which  Is  a large  shallow  estuary  on  the  Texas 
Coast.  The  determination  of  the  minimum  cost  of 
phosphorus  treatment  based  on  source  projections 
for  1990  and  the  predicted  hydrodynamics  (4)  of  the 
bay  are  presented.  The  convergence  of  the  method 
In  the  positive  definite  case  was  demonstrated  bv  a 
simple  ratio  test  for  the  sequence  of  Initial  value 
corrections  generated  by  the  shooting  algorithm. 

The  extension  of  this  approach  shows  promise  In  the 
solution  of  boundary  value  problems  In  which  more 
than  one  Independent  variable  Is  considered. 

Multiple  .Shooting  Solutions  of  a Distributed  System 

The  class  of  problems  discussed  here  have  a state 
equation  given  by  an  elllpUc  partial  differential 
equation  satisfying  mixed  Dlrlchlet  and  Neumann 
boundary  conditions. 


i 

14 

1 

i 


I 

y. 

h 


^ (v(x,y) 


by 


by  y ^ 

C^ -K(x,y)  C = f(x,y,a),  (x,y)«t. 

U) 


subject  to  boundary  conditions 

C(x,y)  - given  (x,y)  t bC,  , Dlrlchlet  condition 

‘ (2) 


Q • n « 0 


(x.yjcbCi^  .Neumann  condition 


(3) 


where  C(  . , . ) - 

E (x,y),E  (x.y)  *(!*'  ^0) 
X y 


f(x,y,-) 

K(x,y) 


state  variable 


u(x,y),  v(x,y) 


(x.y) 


-dispersion  coeffi- 
cients In  the  co- 
ordinate directions, 

E ,E  >0,  V(x,y)  sfl 
X y 

-control  functions 

-reaction  rates 

-outward  pointing 
normal  on  the 
boundary 

-[(uC-E^bC/bx)  , 

(''C-Ey-lf  )] 

-velocity  coefficients 
In  the  coordinate 
directions 


-control  parameters, 
o < R 

-simply,  connected 
open  region , 0 c 

-boundary  of  region  0 


b0»  bfljUbOj 

If  In  addition  to  the  conditions  specified  above,  the 
equations  satisfy  an  elliptlclty  condition  (1),  then 
there  exists  a unique  solution 

C(x,y)  « . 

Furthermore,  the  problem  Is  well-posed  with  respect 
to  the  Initial  data  In  the  sense  of  Hadamaard  (6). 


The  'solution  technique  discussed  In  this  paper  Is 
based  on  reformulating  the  boundary  value  problem 
as  an  Initial  value  problem.  Although  the  resulting 
initial  value  problem  violates  the  Hadamaard  con- 
dition, a special  variant  of  multiple  shooting  Is 
employed  to  Improve  convergence  and  generate  the 
initial  conditions  compatible  with  the  specified 
data  and  the  solution. 

Multiple  Shooting  Solutions 

For  simplicity  consider  the  pure  Dlrlchlet  problem 
defined  on  a convex  region  0.  Define  a partition 
in  the  y-coordlnate  direction  and  approximate  the 
derivatives  along  the  partitioned  coordinate  axis. 

The  approximation,  Cj(x),  then  satisfies  the  two- 

point  boundary  value  problem  defined  for  a sei  of 
coupled  second  order  ordinary  differential  equations 
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subject  to  boundary  corxJltlons 


(5) 


-Independent  variables  where 


X,  , - X coordinate  for  the  beginning  of  the 
Jth  filament 

X,  - X coordinate  for  the  eixl  of  the  jth 
filament 

by  - space  mesh 

Ej^  - E(x,(J+i)by) 

- V(x,  0+1)  by) 
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The  above  equations  (4)  may  be  rewritten  In  a 
blocked  matrix  form  as  follows 


the  deviation;  between  specified  and 
calculated  terminal  conditions  at  the 
boundary  points  given  by 


C (x  ,g  ) - C (x  ,g) 

mAv  * ' 


(6)  The  sequence  of  iterates  is  determined  by  the  formula 

,9) 


[c,(x) c (x)r 

i m 


where  g'  ' Is  arbitrary. 

In  the  Ideal  situation,  then 


[ 6c,(x) ftc  (x)  r 

i m 


and  the  algorithm  converges  In  one  step. 


trldlagonal  matrix  dependent  on  the 
coefficients 


= dlag  [ » , 6, i ];  4 =(u,(x) 


m^'  -J  J' 


I = m x m Identity  matrix 

subject  to  the  boundary  conditions 


^ <’'max)“5‘''®" 

The  problem  is  a linear  two  point  boundary  value 
problem,  which  In  the  Ideal  case  can  be  solved  by 
shooting  methods  (3,7).  The  Instability  of  the 
Initial  value  problem  for  the  elliptic  system  remains 
and  Is  manifested  by  the  rapid  growth  of  the  trial 
solutions  (7)  made  with  Incorrect  Initial  conditions. 

Let  the  missing  Initial  values  be  designated  by 

. I- 1.2.3 m 

Then, the  shooting  solution  corresponds  to  Newton's 
solution  for  a linear  algebraic  system 

? (g)  - ♦ (g  - g ) (8) 

where  ♦ Is  the  Jacobian  matrix,  4,,  “ ^ ^1 


arbitrary  Initial  value  vector,  g c R 

Initial  value  vector  corresponding  to 
the  solution  satisfying  the  boundary 
conditions 


In  practice,  the  shootlrm .algorithm  generates  a 
sequence  of  Iterates  (g  *}  , which  Is  convergent  If 
the  matrix  Y Is  positive  definite  and  It 

1 I ♦ - I I L < 1 (10) 

Due  to  the  growth  of  trial  solutions,  condition  (10) 

may  not  be  satisfied  and  the  ordinary  shooting 

algorithm  would  consequently  fall  to  converge.  The 

effect  of  the  excessive  growth  is  usually  manifested 

by  the  loss  of  significant  figures  in  the  trial  solutions 

but  It  may  be  reduced  by  multiple  shooting  (3,7,8) 

methods.  For  multiple  shooting,  the  strips 

(*-,i„  i)  partitioned  Into  subintervals 

min , j mdx , j 

and  numerous  alternative  multiple  shooting  strategies 
may  be  Implemented.  The  size  of  the  subintervals  Is 
chosen  to  optimize  computational  accuracy  and  en- 
hance the  convergence  of  the  shooting  algorithm. 
Although  the  dimension  of  the  Initial  value  vector,  g. 
Is  increased  by  multiple  shooting,  the  Jacobian  matrix 
Is  very  often  sparse  and  may  be  Inverted  rapidly  by 
special  algorithms.  The  method  will  be  Illustrated 
by  a practical  example,  where  the  initial  value  pro- 
blems are  solved  by  a simple  modified  Euler  numerical 
integrator  such  that  the  resulting  equations  are 
equivalent  to  ordinary  finite  difference  approximation 
of  the  problem.  The  order  of  the  Jacobian  for  the 
multiple  shooting  solution  technique  Is  still  con- 
siderably less  than  the  finite  difference  approximation 
matrix  and  achieves  the  same  order  of  convergence 
accuracy  In  less  than  half  the  computational  time  of 
the  successive  over-relaxation  (4)  method  and  In 
about  one  hundredth  the  number  of  Iterations.  This 
method  has  several  promising  extensions  such  as  the 
variable-metric  multiple  shooting  technique,  which 
eliminates  the  need  to  Invert  the  Jacobian  (9).  The 
applications  of  multiple  shooting  combined  with 
special  fixed  step/ variable  step  numerical  Integrator 
schemes  also  provides  a method  of  decreasing  the 
truncation  errors  in  the  numerical  approximations 
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of  the  state  equations. 


subject  to  boundary  conditions. 


Point  Control  and  Optimization 

An  important  design  problem  Involves  the  selection 
of  constrained  control  parameters  a,  which  achieve 
some  specified  goal  at  minimum  cost,  where  the 
goals  are  assumed  to  lie  within  the  reachable  set 
of  the  admissible  region.  The  problems  of  this 
type  can  be  formulated  as  the  minimization  of  a 
convex  cost  function  subject  to  linear  constraints. 
The  linear  constraints  are  determined  from  the 
sensitivity  functions  to  the  a parameters  and  con- 
stitute part  of  a nonlinear  programming  solution. 
The  large  dimensionality  of  the  constraints  is  re- 
duced by  exploiting  some  of  the  special  qualities 
of  the  sensitivity  functions  which  we  derive  below. 

The  control  problem  may  be  formulated  as  follows 


Sjj(x,y;  = given  , (x,y)  t 60^ 


. n = 0 , (x,y)  t 


where  C^(x,y)  Is  the  solution  to  (1),  (2)  and  (3)  for 
a reference  or  . Assume  that  f Is  linearly  dependent 
on  the  paramiter  vector  a,  then  the  sensitivity  func- 
tion satisfies  a mo.iotonlclty  condition  outlined 
In  Theorem  1 . 

Lemma  (11) : If  f:  D -•  R Is  continuous  and  the  do- 
main D Is  a compact  subset  of  a topological  space, 
then  f attains  Its  maximum  and  Us  minimum,  l.e., 
there  exists  points  Xj , Xj  of  D such  that 

f(Xj)  = Inf  f(x)  and  f(x„)  = sup  f(x) 
xtD  xtD 


Inf  J(a)  : J(a  ) < J(a ),  v ot  e > 
or  c J> 


where  R is  the  set  of  all  real  numbers. 


subject  to  the  constraints  equation  (1),  (2)  and  (3)  Theorem  I: 


J'  c;  R , Is  a convex  set  of  control  parameter 


Consider  the  point  controlled  distributed  parameter 
system  described  by  equations  (1)  and  conditions 
(2)  and  (3)  with  K > 0.  The  sensitivity  functions 


constraints  and  such  that 
C(x,y)  < g(x,y) 

where  for  simplicity  It  Is  assumed  that  the  function 
g(x,y)  c H (Jj,  the  reachable  set  for  the  parameter 
control,  and  the  cost  functional  J(*)  are  assumed 
to  be  convex.  The  problem  defined  above  may  be 
reformulated  as  follows: 


inf  J(a) 
a c i> 

subject  to  constraints. 


(x,y)  + Sj,  • Of  - C(x,y)  < g(s,y) 


depend  monotonlcalty  on  a and 


Sgn  {S^(x)}  “ - Sgn  {— ) . 


provided  that  the  sign  of  r — Is  either  positive  (or 

0 or 


negative  everywhere 


The  sensitivity  function  S (10)  can  be  shown  to 
or 

satisfy  (12)  and  conditions  (13)  and  (14).  Multiply 
(12)  by  s Integrate  over  fl  and  applying  Green's 

Ct 

formula  gives 


^ 6S^  bS„  6S„ 

r IV  o p Of  Of 

JflV^x  6x  6x  * y By  By 


where  S Is  the  sensitivity  function  (10)  matrix 

''be 

(vector)  -r — which  satisfies  the  problem 

BOfj 


B(S  u)  6(S„v)  ^ 


b-tCV’*'’'’  -bT  * T7  6T)‘  ‘“V 

■ -K<x,y)s„»|^(x,y,0f),  (x,y)cn^^^^ 


li;  . (t2'2  (n)  (16) 


From  the  elUptlclty  of  the  system  equations,  the  left 
side  of  (16)  Is  less  than  zero,  hence 


-A 

M 


f 

4.r 

f1 


t 

V 


■ii' 


Consider  the  special  case  when 

= h(x,y)  >0,  V (x,y)  » n , then  equation  (17) 
0 or  — 

Implies ; 


(18) 


In  dxdy  - dxdy  + J dxdy 

p N 


where  0 = {(s,y)cn:S  >0} 

p a — 


( (x,y)  * 0 : S <01 

Of 


It  Is  clear  from  equation  (17)  that  0.,  / i . 

N 

2 

Applying  the  Lemma  to  any  subset  of  0^  c p , It  Is 
clear  that  S attains  a maximum  and  there  exists 

Of 

points 


(Xjj.yg)  € Op  such  that 


bS 



h X 


2 2 
6 S t 

F7'°  and  -^,-^<0 


bS 


(19) 


b X 


by 


If  the  maximum  occurs  on  the  Neumann  boundary, 
then  from  the  condition  (14)  It  follows  that 


bv 


of  Inputs,  and  the  reachable  set  defined  on  the  ad- 
mlssable  region  may  be  readily  expressed  as  the 
range  of  the  affine  mapping  (11a). 

The  result  Is  very  useful  In  the  example  problem  be- 
cause It  provides  a means  of  simplifying  the 
optimization  task  and  It  also  provides  a method  of 
Identifying  the  domain  of  Influence  for  different  con- 
trol sources  In  the  spatial  domain. 

The  approximate  sensitivity  functions  S (1,))  may  be 
determined  directly  by  solving  the  sens^lvlty 
equations  (12)  through  (14)  using  a shooting  algorithm 
similar  to  that  used  for  the  system  state.  An  alterna- 
tive solution  technique  employs  numerical  differentia- 
tion with  respect  to  the  parameters  and  determines 
the  sensitivity  functions  from  a series  of  numerical 
solutions  of  the  original  problem  as  follows: 


Step  I: 
Step  II: 


Solve  the  system  equations  with 
nominal  control  and  determine  C°(x,y). 


Repeat  Step  I,  m times,  each  time  per- 
turbing one  compMent  of  the  control 
parameter  and  generate  the  set  of 
solutions 


where  m Is  the 
c 


C^x,y)  , )-l,2, . . ,m^ 

member  of  controllable  sources.  Choose 
A0(,  the  perturbation  parameter  to  ensure 
that  the  change  In  the  solutions  Is  sig- 
nificantly more  than  the  convergence 
criterion . 


0 , where  v Is  the  coordinate  along  the  nor-  Step  III:  Calculate  the  Influence  matrix  (vector) 


mal.  Applying  the  Lemma  to  the  one  dimensional 
manifold  glv 

Note  that  S 


manifold  given  by  also  gives  (19). 


(x,y)  - C^x,y)  - C°(x,y) 


(21) 


fio 


0 on  bOj  . 


Applying  (19)  to  equations  (12)  and  assuming  that 


b u b V ,,  , 

- — , -7—  are  small  gives 

b X b y 


S < 0 at  (x„,y.)  € f!  , which  contradicts  the 
or  U U p 

assumption  that  (Xg,yg)  t Op  , hence  Op  - 0 . 


Repeating  the  argument  with  the  sign  of  g(x,y)  re- 
versed gives  the  result 


Step  IV: 


The  solution  for  small  changes  (or  In  the 
linear  case)  about  the  nominal  control 
can  be  determined  from 


C(x,y,af)«C(x,y,af.)  + S (x,y)‘(or-o J 
u or  O 


(22) 


The  control  problem  may  then  be  reformulated  as  the 
following  nonlinear  programming  problem  consisting 
of  a (cost)  performance  Index,  linear  Inequality  con- 
straints and  bounds. 


minimize  I(a) 


(23) 


Sgn  { S^(x,y)  1 “ - Sgn  { 1 


(20) 


Q.  E.  D. 


The  above  result  demonstrates  the  monotonlclty  of 
the  sensitivity  function  for  a very  special  class 


subject  to 

where  a c > , the  set  of  admissible  parameters 
and  (x  ,y  ) are  the  coordinates  of  the  finite  difference 
grid  points 
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computer  system.  The  solution  was  achieved  In 
six  Iterations  when  the  convergence  criterion  was 
specified  to  be 

I (n)  (n-1)  I ^ , .-6 

max  1 - g^  1 < 10 

The  convergent  algorithms  were  characterized  by  the 
property  that  the  ratios 

1 9j  9j  1 

were  constant  (<1).  Attempted 


Although  the  number  of  constraints  In  (24)  could  be 
very  large,  a significant  order  reduction  may  be 
achieved  by  using  a set  of  relaxed  constraint 
functions  g^'iXj.y^)  instead  of  g(Xj,yj),  where 

g^(Xj.yj)  tend  to  gCx^.y^). 

The  above  results  are  Illustrated  In  the  following 
water  quality  examples. 

An  Application  of  Multiple  Shooting  Models  In 
Water  Quality  Analysis  and  Management 

The  computation  of  concentrations  for  water 
quality  constituents  In  Corpus  Chrlstl  Bay,  (Figure 
1)  Illustrates  an  Important  application  of  a multiple 
shooting  solution  for  elliptic  distributed  parameter 
systems.  Corpus  Chrlstl  Bay  Is  a large  shallow 
estuary  on  the  Gulf  of  Mexico  coast  of  Texas  with 
an  approximate  area  of  two  hundred  and  sixty  square 
miles.  The  two  dimensional,  time  ayeraged  hydro- 
dynamic  flow  components  u(x,y),  y(x,y)  In  the  bay, 
were  determined  In  an  earlier  study  (12).  The  point 
sources  In  this  problem  are  pollution  discharges 
from  yarlous  Industrial  complexes  located  on  or 
within  the  bay  perimeter.  The  boundary  conditions 
are  specified  by  the  observed(Dlrlchlet)  concen- 
trations at  the  Laguna  Madre,  Aransas  Pass  and 
Rockport  Inlets,  (Fig.  1)  and  the  no  flow  (Neumann) 
conditions  at  the  barriers  and  land/water  Interfaces. 

The  numerical  approximation  technique  utilized  In 
this  example  Is  similar  to  that  discussed  earlier, 
howeyer,  the  numerical  Integration  Is  perfomed  by 
using  a second  order  approximation  for  equations 
(4),  The  resulting  system  Is  equivalent  to  the 
finite  difference  approximation  defined  on  the 
grid  shown  In  Figure  2,  Numerous  alternative 
shooting  policies  could  be  defined  for  this  pro- 
blem. The  directions  of  the  shooting  computations 
used  In  the  example  are  Indicated  In  the  diagram 
(Fig.  2)  . For  the  case  Illustrated,  32  Initial  and 
Intermediate  values,  indicated  as  on  Fig.  2, 

were  necessary  to  completely  reformulate  the 
resulting  system  as  an  explicit  initial  value  pro- 
blem. A consistent  set  of  terminal  conditions  In- 
dicated as  T's  on  Figure  2 were  determined  by 
checking  the  continuity  and  consistency  of  the 
solutions  when  computed  from  different  directions 
as  Indicated  In  Figure  2. 

The  system  was  solved  for  the  phosphorus  con- 
centrations throughout  the  bay.  The  measured 
and  calculated  values  at  16  observation  points 
throughout  the  bay  were  compared.  The 
phosphorus  distribution  corresponding  to  the 
minimum  square  error  between  observed  and  cal- 
culated values  (Table  1)  is  shown  In  Fig.  3.  The 
solution  for  the  multiple  shooting  case  was  com- 
puted In  less  than  6 seconds  on  a CDC6600 


acceleration  of  the  solution  by  extrapolation  methods 
using  the  above  ratios  slowed  convergence  because 
of  the  effect  of  round  off  on  the  computations.  The 
errors  (Figure  4)  can  be  seen  ^o  be  reduced  Initially 
by  a very  large  value  of  ^10  and  then  by  smaller 
factors  and  finally  oscillations  due  to  roundoff.  The 
number  of  Iterations  required  for  similar  accuracy 
for  a method  equivalent  to  the  successive  over- 
relaxation  technique  (4)  was  about  five  hundred  and 
required  about  14  seconds  on  the  CDC6600  compu- 
ter. It  Is  clear  from  this  study  that  multiple 
shooting  methods  are  competitive  with  and  in  some 
cases  better  than  other  well  known  Iterative  methods. 
The  advantages  Include  the  reduced  dimension  of  the 
matrices  to  be  solved,  which  In  this  example  were 
solved  by  a standard  linear  equation  subroutine. 

For  larger  problems,  the  solutions  may  be  computed 
by  combining  a variable  metric  technique  with  the 
shooting  method  to  eliminate  some  of  the  Inversions 
or  linear  system  solutions. 


Determine  the  optimal  phosphorus  treatment  policies 
at  the  sources  for  some  target  year  In  the  future 
such  that  the  total  cost  of  keeping  the  water  quality 
at  1969  levels  Is  minimized.  This  problem  may  be 
equivalently  stated  as  follows 
M 

min  I(or)  “ Coat  (Oj)  (25) 


subject  to 

Sj,(lJ)  «<  Cjggg  -Cg(l.j)  (26) 

where  the  cost  function  for  a chemical  coagulation 
plant  assuming  a twenty  year  life  Is 


The  resulting  model  can  be  used  to  evaluate  treat- 
ment policies  for  the  discharges  Into  the  estuary, 
and  a typical  water  quality  management  problem 
may  be  stated  as  follows; 


Cost  (Oj)  -SxlO^k^O  ^ + 0.365xlo\^Q  ^ $/year(4) 

(27) 


0 - , 


$264, 000/year  are  as  shown  In  Table  1 with  the 
corresponding  optimal  profile  shown  In  Figure  7. 
This  corresponds  to  an  overall  savings  of  about 
ninety  thousand  dollars  per  year  over  the  cost  of 
maintaining  total  effluent  discharges  for  the  eleven 
controllable  sources  at  the  1969  levels. 

SUMMARY  AND  CONCLUSIONS: 


M = number  of  controllable  sources  (11) 

Tlj  = 0.93:  = 0.098:  = 0.90 


k,  = 5.45:  k.  = -0.03. 

3 4 


Qj  = size  of  the  controllable  discharge 


where  the  Indices  1 and  ) correspond  to  the  grid 
point  source,  and 


[a  , J=1 , 2, . , . ,m  : 0,01<o  <0.93} 
J c j — 


The  parameter  sensitivity  function  matrix  for  1960 
conditions  Incorporating  source  (4)  and  hydro- 
dynamic  projections  for  Corpus  Chrlstl  Bay  was 
determined  by  the  numerical  differencing  algorithm. 
The  sensitivity  function  to  the  source  Intensity 
for  one  of  the  pollution  discharges  Is  shown  In 
Figure  5 and  the  percent  difference  between  the 
concentrations  predicted  by  the  shooting  models 
and- the  sensitivity  model  Is  shown  In  Figure  6. 

The  approximate  sensitivity  functions  In  general 
satisfy  the  condition  specified  In  equation  (17), 
however,  some  of  the  functions  were  positive 
in  the  upper  regions  of  the  bay  (Nueces  Bay  and 
River),  where  the  positive  definiteness  of  the 
operator  was  questionable.  The  sensitivity  model 
Is  within  1%  of  the  shooting  model  except  In  the 
small  area  representing  Corpus  Chrlstl  harbour 
where  the  convergence  error  for  the  shooting 
algorithm  was  larger  than  usual.  The  large  error 
observed  In  the  upper  regions  of  the  bay  Is  due 
to  loss  of  numerical  accuracy  due  to  the  growth  of 
the  Initial  value  solutions.  This  type  of  error  can 
be  eliminated  by  a shooting  scheme  Involving 
shooting  In  the  negative  y direction  starting  at  the 
Nueces/Corpus  Chrlstl  Bay  Junction. 


A method  of  solving  the  sparse  large  linear 
approximations  for  point  controlled  elliptic  dis- 
tributed parameter  systems  by  multiple  shooting  was 
discussed.  The  method  was  applied  to  the  analysis 
of  a water  quality  system  with  264  elements.  The 
linear  system  generated  by  the  multiple  shooting 
technique  was  a 32x32  matrix  with  a large  number 
of  zero  elements.  Convergence  was  achieved  In  a 
few  Iterations  at  a speed  twice  that  of  a compara- 
tive successive  over-relaxation  algorithm  for  the 
same  problem.  A result  about  the  sensitivity  func- 
tion of  the  control  parameters  was  derived,  which 
serves  as  a useful  check  for  an  alternative  per- 
turbation model  for  the  problem.  An  analysis 
technique  based  on  the  shooting  algorithm  was  em- 
ployed to  solve  for  the  optimal  management  policies 
for  controllable  sources  In  Corpus  Chrlstl  Bay.  The 
method  may  be  extended  to  solve  even  larger  pro- 
blems by  a combination  of  multiple  shooting  with 
the  variable  metric  methods  which  eliminate  the 
need  to  solve  large  linear  algebraic  systems. 
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TABLE  1 

Observed  and  Calculated  Long-Term 
Phosphorous  Concentrations  In  Corpus  Chrlstl  Bay 


1 j 


TABLE  2 

Summary  of  Optimal  Treatment 


Controllable 
Source  No. 

Plant  Location 
Indicles  (1,J) 

Plant  Flow 
(MGD) 

Percent  Removal 
Efficiency 

Treatment  Cost  In 
Dollars 

(Average,  per  year) 

3 

15,3 

2.779 

4.18 

21,927 

6 

20,7 

0.646 

17.5 

21,926 

7 

28,7 

0.582 

18.7 

21,926 

10 

4,9 

10.340 

6.1 

23,586 

11 

2,10 

4.265 

4.0 

21,996 

12 

14,10 

0.194 

67.4 

21,937 

13 

13,  14 

1.680 

6.8 

21,926 

14 

7,16 

8.725 

62.9 

43,058 

16 

12,16 

0.969 

14.2 

21,944 

17 

7,17 

131.129 

0.1 

21,932 

23 

9,20  1 

4.265 

2.9 

21,926 

Total  Cost  = 264, 0B4  Dollars/year 
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